If anything remains unclear or unanswered after this session, also if you spot mistakes or typos, contact:
Twitter: @liza_p_semenova
Github: elizavetasemenova
E-mail: elizaveta.p.semenova@gmail.com
Overview:¶
Julia:
- Why Julia?
- basic syntaxis
- How to continue learning Julia after this workshop
Bayesian Inference:
- What is Bayesian inference?
- Why do we care?
- How to perform Bayesian inference: a painful way and a modern way
- Examples
Outro
- Bayesian Neural Networks
- Bayesian decision theory
- Closing remarks
We will use Julia for the practical part:
println("All models are wrong, but some are useful.")
Docker¶
Follow instructions here: https://github.com/elizavetasemenova/EmbracingUncertainty
Also Python and R versions of the Bayesian block are available in this repository.
Julia¶
Why Julia¶
- Pyhton and R are easy to learn but lack speed on native code.
- Julia was created (released in 2012) to overcome this issue.
- Julia is not meant to substitute the two established GPPLs (general purpose programming languages), but rather to augment them since a lot of R and Python functionality is available in Julia.
- Julia's popularity is growing in the scientific community:
https://www.nature.com/articles/d41586-019-02310-3
- Julia provides a fertile ground for developement of probabilistic programming languages (PPLs): several native ones have been recently developed (Turing, Gen, Stheno, SOSS, Omega) indicating strong interest of the Bayesian research community.
Basic Julia syntaxis¶
##############################################
# basic Julia syntaxix
##############################################
println("Hello, world!")
3 - 2.4
println("computed ", 3 - 4.5)
println("Hello" * " Julia!")
" who" ^ 3 * ", who let the dogs out" # do you know this song?
Arithmetic operations¶
- x + y binary plus
- x - y binary minus
- x * y times
- x / y divide
- x ÷ y integer divide
- x \ y inverse divide
- x ^ y power
- x % y remainder
println((50 / 60) ^ 2)
x = 9
y = 7
x % y
rem(9, 7)
x ÷ y
x \ y
Logical expressions¶
- ! is the NOT operator
- && is the AND operator
- || is the OR operator
Beware:
- & is the "bitwise AND"
- | is "bitwise OR"
false && true || true
Julia allows one to omit the multiplication operator when writing expressions.
x = 3
4x^2+x
4+1<2||2+2<1
a=true
b=!!!!a
Greek letters :¶
type in '\theta'
This allows to copy math equatios directly, as if we would be writing them by hand.
θ = 3
2θ
Type system¶
typeof(1)
typeof("Hello, world!")
typeof("αβγδ")
typeof("H")
typeof('H')
typeof(1+3+5.)
Array types¶
Array{Int64}(undef, 3)
Array{String}(undef, 3)
'a'==="a"
x=3
y=2
a = Array{Integer,2}(undef, x, y)
x = Array{Int64}(undef,11, 12)
typeof(x)
Functions¶
a, b, c = cos(0.2), log(10), abs(-1.22) # multiple assignment
myfunc(x) = 20*x
myfunc(2)
add(x,y) = x + y
add(33, -22.2)
function nextfunc(a, b, c)
a*b + c
end
nextfunc(7,5,3)
function print_type(x)
println("The type of testvar is $(typeof(x)) and the value of testvar is $x")
end
a = ['1',2.]
print_type(a)
function f(x)
return 2x
3x
end
f(5)
Multiple methods¶
cos_func(x) = cos(x)
cos_func(.7)
cos_func(adj, hyp) = adj/hyp
cos_func(.7)
cos_func(12, 13)
methods(cos_func)
?cos_func
cos_func(theta::Float64) = cos(theta) # :: forces Julia to check the type
Use '?' to searach. For example, '?cos':
?cos
Working with arrays¶
x = [1, 6, 2, 4, 7, 2, 76, 5]
size(x)
length(x)
typeof(x)
x = collect(1:7)
x = range(0, stop = 1, length = 10)
x = collect(x)
x = fill(5, 4)
Vectorized "dot" operators¶
x + 2
x .+ 2
x .* x
data = [1.6800483 -1.641695388;
0.501309281 -0.977697538;
1.528012113 0.52771122;
1.70012253 1.711524991;
1.992493625 1.891000015]
size(data)
typeof(data)
data[:,1] # all rows, 1st column
data = [[3,2,1] [3,2,1] [3,2,1] [3,2,1] [3,2,1] [3,2,1] [3,2,1] [3,2,1] [6,5,4]]
data[2,:]
rows, cols = size(data)
println(rows)
println(cols)
Loops¶
for num = 1:2:length(x)
println("num is now $num")
end
values = [ "first element", 'θ' , 42] # an array
for x in values
println("The value of x is now $x")
end
x = [3, 2, 1]
count=1
for i in x
x[i]=count
count=count+1
end
println(x[3])
println(x)
if 5>3
println("test passed")
end
Plots¶
using Plots
gr()
x = cos.( - 10:0.1:10)
plot(1:length(x), x)
plot!(title = "First plot", size = [300, 300], legend = false) # modify existing plot
hline!([0])
plot(1:length(x), x,
size = [300, 300],
legend=false,
line= :scatter,
marker = :hex)
x = collect(1:0.1:7)
f(x) = 2 - 2x + x^2/4 + sin(2x)
plot(x, f)
Save plot¶
#savefig("Plot_name.png") # Saved png format
Distributions.jl¶
using Distributions
d = Normal()
n_draws = 1000
draws = rand(d, n_draws)
histogram(draws, size = [300, 300], bins = 40, leg = false, norm = true)
[pdf(d, x) for x in -3:1:3] # evaluate PDF/PMF
shape = (3, 2)
rand(Uniform(0,1), shape)
How to continue learning Julia after this course¶
Please run these commands now:
using Distributions
using Plots
using StatsBase
using Turing
Bayesian Inference¶
WHAT¶
What is Bayesian inference?¶
Bayesian inference is the process of deducing properties of a probability distribution from data using Bayes’ theorem.
Bayes' theorem relies heavily on the notion of conditional probability.
Conditional probability¶
Conditional probability of A given B:
$$P(A \vert B)=\frac{P(A \cap B)}{P(B)}.$$
$P(A \vert B)$ - probability of event A given B has occured
$P(A \cap B)$ - probability of event A occured and B occured
We can derive
The conditional probability of A given B is the conditional probability of B given A scaled by the relative probability of A compared to B.
Bayes' Theorem example 1: diagnosing flu¶
Find probability of a patient having a flue given they have high temperature.
A = to have a flu
B = to have high fever
$$P(\text{flu} \vert \text{fever})= P(\text{fever} \vert \text{flu}) \frac{P(\text{flu})}{P(\text{fever})}$$We know that the probability of having fever this time of the year is 10%, probability of having a flue this time of the year is 7%, and among all people having a flu, 70% of them have fever.
$$ 70\% * 7\% / 10\% = 49\%.$$Bayes' Theorem example 2: spam filtering¶
We would like to test a message for being a spam, by checking that it
contains a predefined set of words or word combinations, which are
common for spam messages ("Free investment", "Dear email"):
$$P(\text{spam } \vert \text{ words}) = \frac{ P(\text{words } \vert \text{ spam}) P(\text{spam})}{P(\text{words})}$$
In practice, we usually want to estimate the unknown model parameters $$\theta = (\theta_1, ..., \theta_k),$$ given the observed data $$y = (y_1, ..., y_n)$$ by using the Bayes' rule: $$P(\theta \vert y) = \frac{P(y \vert \theta) P(\theta)}{P(y)}.$$
In continuous case, we speak of probability densities $f(.)$.
The posterior distribution is being formed as
$$f(\theta \vert y) \propto f(y \vert \theta) f(\theta).$$
$f(\theta)$ is the prioir and expresses our beliefs about $\theta$
$f(y \vert \theta)$ is the likelihood of data $y$ as a function of $\theta$
we dropped $f(y)$ because the average likelihood across all parameter values $$f(y) = \int_{\theta_1} ... \int_{\theta_k} f(\theta) f(y \vert \theta) d\theta_1 ... d\theta_k = E_{\theta}f(y \vert \theta)$$ does not depend on $\theta$ (hence will not influence the inference) and is very hard to compute.
Components of the bayesian model¶
Likelihood
Prior
Posterior
WHY¶
Diagnosing the cause of headache¶
Frequentist doctor:
- A mental model for the cause of pain
- Perform tests
Bayesian doctor:
- A mental model for the cause of pain
- History of the patient
- Perform tests
I.e. frequentist inference is like a doctor who never reads patient's history even for patients with chronic conditions.
Frequetist vs Bayesian approaches¶
Frequentist:
- probability is a long-run frequency
- parameters are fixed,
- parameters can only be estimated using the current data.
Bayesian:
- probability represents one's uncertainty
- parameters are described by distributions (i.e. reflecting uncertainty),
- estimation is based on our belief and/or prior knowledge about parameters and the current data.
Remark: historically, Bayesian approach appeared earlier than the frequentist one.
Why have we been using frequentist methods all these years?
Because we didn't have computers! Bayesian inference is very hard to make by hand.
Why use Bayesian inference?¶
Prior knowledge
Flexible in uncertainty modeling, particularly under small amount of available data
Very flexible model formulation accounting for the mechanistic knowledge about a system
Note:
- Results are essentially equal with enough data, i.e. the posterior approximates the MLE
Exercise¶
- think of a problem relevant to your work where one could benefit from the Bayesian approach
HOW¶
How¶
analytically
elegant, but rarely possible.
numerically
Create posterior samples, describing the distributions of parameters.
We achieve this by exploring the space of parameters to find the most probable combinations of parameters.
Further we treat the obtained sampled as new data, to extract information about parameters, such as mean, credible interval or other statistics.
Analytically¶
Numerically¶
- Markov Chain Monte Carlo (MCMC) family of algorithm, e.g.
- Metropolis-Hastings
- Gibbs
- Hamiltonian Monte Carlo (HMC)
- No-U-Turn sampler (NUTS)
- further variants such as SGHMC, LDHMC, etc
- Variational Bayes
- Approximate Bayesian Computation (ABC)
- Particle filters
- Laplace approximation
Priors¶
informative / non-informative
conjugate priors guarantee that the posterior has an easily calculable form
Some conjugate families of distributions
(Robert, Casella, "Monate Carlo Statistical Methods")
Example - Priors¶
using Distributions
using Plots
using StatsBase
##############################################
# prioir x likelihood = posterior
##############################################
success=6
tosses=9
# Create a distribution with n = 9 (e.g. tosses) and p = 0.5.
d = Binomial(tosses, 0.5)
pdf(d, success)
# define grid
grid_points = 100
p_grid = range(0, stop = 1, length = grid_points)
# compute likelihood at each point in the grid
likelihood = [pdf(Binomial(tosses, p), success) for p in p_grid]
# define prior
prior = ones(length(p_grid));
# As Uniform prior has been used, unstandardized posterior is equal to likelihood
# compute product of likelihood and prior
posterior = likelihood .* prior;
function computePosterior(likelihood, prior)
# compute product of likelihood and prior
unstd_posterior = likelihood .* prior
# standardize posterior
posterior = unstd_posterior / sum(unstd_posterior)
p1 = plot(p_grid, prior, title = "Prior")
p2 = plot(p_grid, likelihood , title = "Likelihood")
p3 = plot(p_grid, posterior, title = "Posterior")
plot(p1, p2, p3, layout=(1, 3), label="")
end
prior1 = ones(length(p_grid))
posterior1 = computePosterior(likelihood, prior1)
Posterior is the same as the prior.
#prior2 = 2 * (p_grid .>= 0.5)
prior2 = 0.5 * (p_grid .>= 0.5)
posterior2 = computePosterior(likelihood, prior2)
Posterior vanishes at the points where prior vanishes. I.e. we need to be careful and avoid assignign zero prior values to the regions where have no evidnece of those values being impossible.
prior3 = exp.(-5 * abs.(p_grid .- 0.5))
posterior3 = computePosterior(likelihood, prior3)
Some priors, such as exponential and Laplace, allow to obtain the "shrinkage" effect.
The Monte Carlo method¶
is a group of algorithms which use random sampling repeatedly to make numerical estimations of unknown qunatities/parameters.
(originated within the Manhattan Project thanks to Stanislaw Ulam)
Find approximate value of $\pi.$
(Note: we are solving a deterministic problem using random number generation)
##############################################
# the Monte Carlo method - compute pi
##############################################
function in_circle(x, y, r)
sqrt(x^2 + y^2) <= r
end
function approx_pi(r, n)
xs, ys, cols = [], [], []
count = 0
for i in range(1, step=1, stop=n)
x = rand(Uniform(0,1))
y = rand(Uniform(0,1))
append!(xs, x)
append!(ys, y)
if in_circle(x, y, r)
count += 1
cols = vcat(cols, :red)
else
cols = vcat(cols, :steelblue)
end
end
pi_appr = round(4 * count/n, digits = 3)
pl = scatter(xs,
ys,
color=cols,
size=(200,200),
legend = false,
xticks = false,
yticks = false,
framestyle = :box,
title = "pi (approximately) = " * string(pi_appr),
titlefontsize=font(7, "Calibri"))
display(pl)
end
r = 1
n = 100
for n in 5 * 10 .^[1, 2, 3]
approx_pi(r, n)
end
Monte Carlo integration¶
(Again, let's solve a deterministic problem using random sampling)
Find value of the integral
$$\int_a^b f(x)dx. $$Monte Carlo integration estimates this integral by estimaing the fraction of random points that fall below $f(x)$.
In our context, we are interested in estimating expectations
$$ E[h(x)] = \int h(x)f(x)dx,$$which can be done with
$$ \bar{h}_n = \frac{1}{n} \sum_i^n h(x_i),$$where $x_i ∼ f$ is a draw from the density $f$.
The convergence of Monte Carlo integration is $0(n^{1/2})$ and independent of the dimensionality. Hence Monte Carlo integration generally beats numerical intergration for moderate- and high-dimensional integration since numerical integration (quadrature) converges as $0(n^d)$.
Example¶
Estiamte the integral $\int_0^1 e^x dx$.
##############################################
# the Monte Carlo method - integration
##############################################
exp(1) - exp(0)
x = range(0, stop = 1, length = 100)
plot(x, exp.(x), size= [200,200], legend= false)
pts = rand(Uniform(0,1), (100, 2)) # sample uniformly in the square
pts[:, 2] *= exp(1)
cols = fill(:steelblue, 100)
for i in range(1, step=1, stop=100)
if pts[i,2] > exp(pts[i,1]) # acceptance / rejection step
cols[i] = :red
end
end
scatter!(pts[:, 1], pts[:, 2], color = cols, size=[250, 250], legend = false, xlim = [0,1], ylim = [0, exp(1)])
# Monte Carlo approximation
for n in 10 .^[1, 2, 3, 4, 5, 6, 7, 8]
pts = rand(Uniform(0,1), (n, 2))
pts[:, 2] *= exp(1)
count = sum(pts[:, 2] .< exp.(pts[:, 1]))
volume = exp(1) * 1 # volume of region
sol = (volume * count)/n
println(sol)
end
Mandatory coin tossing example¶
Is the coin biased if we observed (H, H, T, H, ... , T, H)?
The frequentist way - analytically¶
How to derive the MLE estimate?
and so
$$ \hat{\theta} = \frac{\text{heads}}{\text{heads + tails}} $$The frequentist way - numerically¶
##############################################
# coin tossing
##############################################
n = 4
h = 3
p = h/n
The Bayesian way¶
Beta distribution:
$$ \text{Beta}_\theta(a,b) = C * \theta^{(a-1)} (1 - \theta)^{(b-1)} $$If we compute the posterior distribution analytically, we will obtain a Beta distribution with new parameters.
a, b = 10, 10 # hyperparameters
prior = Beta(a, b) # prior
post = Beta(h+a, n-h+b) # posterior
function beta_binomial(n, h, a, b)
# frequentist
p = h/n
mu = mean(Binomial(n, p))
# Bayesian
thetas = range(0, stop=1, length=200)
prior = pdf.(Beta(a, b), thetas)
post = pdf.(Beta(h+a, n-h+b), thetas)
likelihood = n * [pdf(Binomial(n, p), h) for p in thetas];
plot(thetas,
prior,
size= [400, 400],
label = "Prior",
color = :blue,
xlim = [0, 1],
xlabel = "theta",
ylabel = "Density")
plot!(thetas, post, label = "Posterior", color = :red)
plot!(thetas, likelihood, label="Likelihood", color = :green, legend = :topleft)
vline!([(h+a-1)/(n+a+b-2)], color = :red, linestyle = :dash, label="MAP")
vline!([mu / n], color = :green, linestyle = :dash, label="MLE")
end
beta_binomial(100, 80, 10, 10)
beta_binomial(4, 3, 10, 10)
beta_binomial(4, 3, 2, 2)
beta_binomial(4, 3, 1, 1)
Markov Chain Monte Carlo¶
- We want to estiamte the posterior distribution, but this is often intractable
Markov Chain Monte Carlo idea:
- draw samples from a (simple) proposal distribution so that each draw depends only on the state of the previous draw (i.e. the samples form a Markov chain)
under certain conditions, the Markov chain will have a unique stationary distribution
we set up an acceptance criteria for each draw based on comparing successive states with respect to a target distribution that enusre that the stationary distribution is the posterior distribution we are searching for
there is no need to evaluate the potentially intractable marginal likelihood
after sufficient number of iterations, the Markov chain of accepted draws will converge to the staionary distribution, and we can use those samples as (correlated) draws from the posterior distribution, and find functions of the posterior distribution
Metropolis-Hastings random walk algorithm¶
Start with an initial guess for $\theta$
Chose a new proposed value as $\theta_p = \theta + \delta_\theta, \delta_\theta \sim N(0, \sigma).$
Here we have chosen the proposal distribution to be $N(0, \sigma).$
If $g$ is the posterior probability, calculate the ratio $\rho = \frac{g(\theta_p \mid X)}{g(\theta \mid X)}$
(adjust for symmetry of the proposal distribution)
- If $\rho \ge 1,$ accept $\theta = \theta_p;$ if $\rho < 1,$ accept $\theta = \theta_p$ with probability $p,$ otherwise keep $\theta = \theta.$ (This step is done with the help of the standard Uniform distribution)
##############################################
# Metropolis-Hastings
##############################################
function target(likelihood, prior, n, h, theta)
if (theta < 0 || theta > 1)
return 0
else
return (pdf(likelihood(n, theta), h) * pdf(prior, theta))
end
end
n = 100
h = 61
a = 10
b = 10
likelihood = Binomial
prior = Beta(a, b)
sigma = 0.3
naccept = 0
theta = 0.1
niters = 10000
samples = zeros(niters+1)
samples[1] = theta
for i=1:niters
theta_p = theta + rand(Normal(0, sigma))
rho = min(1, target(likelihood, prior, n, h, theta_p)/target(likelihood, prior, n, h, theta ))
u = rand(Uniform(0,1))
if u < rho
naccept += 1
theta = theta_p
end
samples[i+1] = theta
end
println("Portion of accepted steps = " * string(naccept/niters))
nmcmc = Int(round(length(samples)/2))
post = Beta(h+a, n-h+b)
thetas = range(0, stop=1, length=200)
histogram(samples[nmcmc:length(samples)] ,
size = [500, 300],
label="Distribution of posterior samples", alpha = 0.5,
legend = :topleft)
histogram!(rand(prior, nmcmc),
label = "Distribution of prior samples", alpha = 0.5)
plot!(thetas, 50*[pdf(post, theta) for theta in thetas], color = :red, label = "True posterior")
We run the chain for $N$ iterations and discard the first $B$ samples. This is called burn-in.
We can run several parallel versions of the algorithm. Each of them is called a chain.
Neigbouring samples will contain similar information. We might want to save only every second, or fifth, or tenth. This is called thinnning.
Convergence diagnostics¶
Rigorous way of assesing convergence is an unsolved problems. But there are several tool swe can use to convice ourselves that an MCMC has converged, such as
- trace plots need to look stationary
- parallel chains should carry similar information
function mh_coin(niters, n, h, theta, likelihood, prior, sigma)
samples = [theta]
while length(samples) < niters
theta_p = theta + rand(Normal(0, sigma))
rho = min(1, target(likelihood, prior, n, h, theta_p)/target(likelihood, prior, n, h, theta ))
u = rand(Uniform(0,1))
if u < rho
theta = theta_p
end
append!(samples, theta)
end
return samples
end
n = 100
h = 61
lik = Binomial
prior = Beta(a, b)
sigma = 0.05
niters = 100
chains = [mh_coin(niters, n, h, theta, lik, prior, sigma) for theta in range(0.1, stop=1, step=0.2)];
Compare multiple chains
p = plot(chains[1], size= [500, 500], legend =:false, xlim = [0, niters], ylim = [0, 1])
for i in 2:length(chains)
plot!(chains[i])
end
display(p)
Was this painful to write a sampler by hand?
If not, bare in mind that we only wrote the simplest one possible! Sampling algorithms can get very complicated. for chain in chains pl = plot(chain) display(pl) end
Probabilistic programming languages (PPLs)¶
Luckily, we do not need to write a sampler by hand every time, because probabilistic programming languages (PPLs) are there to help.
A PPL allows to formalize a Bayesian model and perform inference with the help of powerful algorithms. A user needs to only formulate the model (and maybe chose a sampler) and press the inference button.
The list of currently existing PPLs is overwhelmingly long and only keeps growing:
- BUGS, WinBUGS, JAGS,
- Stan,
- PyMC3, PyMC4,
- Nimble,
- Pyro,
- Edward, TensorFlow Probability, Edward 2,
- Gen,
- Turing,
- Stheno,
- SOSS,
- Omega,
- Infer.NET
to name a few.size(chains)
There is already a couple PPLs accessible via or native to Julia.
We will work with Turing.
Turing¶
##############################################
# Turing
##############################################
using Turing
Turing syntaxis¶
Model creation
@model mod(y) = begin
# model definition
end
Coin tossing problem - Turing¶
n = 100 # number of trials
h = 61 # number of successes
niter = 10000
@model coin(n, h) = begin
# prior
p ~ Beta(2, 2)
# likelihood
h ~ Binomial(n, p)
end
Sampling
ch = sample(coin(n,h), NUTS(niter, 0.65));
show(ch)
# read samples into array
p = convert(Array{Float64}, ch[:p].value.data[:,:,1][:,1]);
histogram(p, size = [300, 300], legend = false, title = "posterior density")
# traceplot
plot(p, size = [300, 300], legend = false, title = "traceplot")
function plot_par(par)
p1 = histogram(par, size = [400, 300], legend = false, title = "posterior density")
p2 = plot(par, title = "traceplot")
plot(p1, p2, layout=(1, 2), label="")
end
plot_par(p)
Hierarchical models¶
Hierarchical models put priors on the hyperparamaters.
$$ y \sim f(y \mid \theta) $$$$ \theta \sim h(\theta \mid \sigma) $$More levels of hiearchy are possible.
In this way a hierarchical model allows information to be shared between parameters $\theta,$ since they might be not independent, but rather drawn from a common distribution.
##############################################
# hierarchical models
##############################################
@model coin_hier(n, h) = begin
# hyperparameters
alpha_hyp ~ InverseGamma(10, 2)
beta_hyp ~ InverseGamma(10, 2)
# prior
p ~ Beta(alpha_hyp, beta_hyp)
# likelihood
h ~ Binomial(n, p)
end
niter = 20000
ch = sample(coin_hier(n,h), NUTS(niter, 0.30));
show(ch)
# read samples into array
p = convert(Array{Float64}, ch[:p].value.data[:,:,1][:,1]);
plot_par(p)
alpha_hyp = convert(Array{Float64}, ch[:alpha_hyp].value.data[:,:,1][:,1]);
plot_par(alpha_hyp)
beta_hyp = convert(Array{Float64}, ch[:beta_hyp].value.data[:,:,1][:,1]);
plot_par(beta_hyp)
Estimating parameters of the normal distribution¶
Data generation¶
##############################################
# normal distribution
##############################################
N = 2000
y = rand(Normal(0,1), N)
histogram(y, size = [300, 300], legend = false)
Turing Model: unknown mean, known variance¶
@model norm_mu(y) = begin
sigma = 1
# prior
mu ~ Normal(0,0.5)
# likelihood
for i in eachindex(y)
y[i] ~ Normal(mu, sigma)
end
end
ch = sample(norm_mu(y), NUTS(niter, 0.65));
mu = ch[:mu].value.data[:,:,1]
plot_par(mu)
Multiple chains¶
num_chains = 4
chains = mapreduce(c -> sample(norm_mu(y), NUTS(niter, 0.65)), chainscat, 1:num_chains)
mu = chains[:mu].value.data
plot(mu[:,:,1], label ="chain 1")
plot!(mu[:,:,2], label ="chain 2")
plot!(mu[:,:,3], label ="chain 3")
plot!(mu[:,:,4], label ="chain 4")
histogram(mu[:,:,1], alpha = 0.5, label = "chain 1")
histogram!(mu[:,:,2], alpha = 0.5, label = "chain 2")
histogram!(mu[:,:,3], alpha = 0.5, label = "chain 3")
histogram!(mu[:,:,4], alpha = 0.5, label = "chain 4")
Turing Model: unknown mean, unknown variance¶
@model norm_mu_sigma(y) = begin
# priors
mu ~ Normal(0,0.5)
sigma ~ InverseGamma(2, 3)
# likelihood
for i in eachindex(y)
y[i] ~ Normal(mu, sigma)
end
end
ch = sample(norm_mu_sigma(y), NUTS(niter, 0.65));
show(ch)
mu = ch[:mu].value.data[:,:,1]
sigma = ch[:sigma].value.data[:,:,1]
pl_mu = plot_par(mu)
vline!([0], color = :green, label="MLE")
pl_sigma = plot_par(sigma)
vline!([1], color = :green, label="MLE")
display(pl_mu)
display(pl_sigma)
Linear regression¶
Likelihood:
$$ y \sim N(ax+b, \sigma^2)$$Priors:
$$ a \sim N(0,100) $$$$ b \sim N(0,100) $$$$ \sigma \sim U(0,20) $$Data generation¶
##############################################
# linear regression
##############################################
n = 100
a_true = 6
b_true = 2
x = range(0, stop=1, length = n)
x = convert(Array, x)
y = a_true*x .+ b_true + rand(Normal(0,1), n);
plot(x, a_true*x .+ b_true, legend = false, size = [350, 350], color = :blue)
scatter!(x, y)
@model lin_reg(x, y) = begin
a ~ Normal(0, 10)
b ~ Normal(0, 10)
lp = a * x .+ b
s ~ InverseGamma(2, 3)
for i in eachindex(y)
y[i] ~ Normal(lp[i], sqrt(s))
end
end
niter = 20000
ch = sample(lin_reg(x, y), NUTS(niter, 0.65));
a = ch[:a].value.data[:,:,1]
b = ch[:b].value.data[:,:,1]
s = ch[:s].value.data[:,:,1]
pl_a = plot_par(a)
vline!([a_true])
pl_b = plot_par(b)
vline!([b_true])
pl_s = plot_par(s)
vline!([1])
display(pl_a);
display(pl_b)
display(pl_s)
Binomial likelihood¶
##############################################
# binomial likelihood
##############################################
function invlogit(x)
exp.(x) ./ (1 .+ exp.(x))
end
n = 1000
x = rand(Normal(), n)
alpha_true = -0.3
beta_true = 0.7
ps = alpha_true .+ beta_true*x
ps = invlogit(ps)
y = [rand(Binomial(10, p)) for p in ps];
histogram(ps, bins = 20, title = "p", label = "", size = [300, 300])
@model binom(x, y) = begin
alpha ~ Normal(0, 1)
beta ~ Normal(0, 1)
p = invlogit(alpha .+ beta*x)
for i in eachindex(y)
y[i] ~ Binomial(10, p[i])
end
end
ch = sample(binom(x, y), NUTS(niter, 0.65));
show(ch)
alpha = ch[:alpha].value.data[:,:,1]
beta = ch[:beta].value.data[:,:,1]
pl_a = plot_par(alpha)
pl_b = plot_par(beta)
display(pl_a);
display(pl_b);
OUTRO¶
Bayesian Neural Networks¶
- Neural networks (the tool of deep learning) is composition of functions on matrices
- Can we supply predictions of neural networks with uncertainty?
- Mind the gap: interpretability of priors
(Graphics by Eric Ma)
- Bayesian deep learning: learning a probability distribution of each parameter, rather than point estimates
(Graphics by Eric Ma)
"Gaussian Process Behaviour in Wide Deep Neural Networks"
Alexander Matthews
R. Neal (1994):
A single layered network (MLP) with $D$ input features and $K$ hidden units, where prior is appropriately scaled, i.e.
- normal distr with 0 mean and variance is $\frac{1}{D}$ for the "input -> Layer 1" weights
- normal distr with 0 mean and variance is $\frac{1}{K}$ for the "Layer 1 -> output" weights
As $K$ tends to infinity, the hidden layer converges in distribution to a multivariate normal $N(0,C)$ with covariance matrix $C.$
(Proof: Mutlivariate version of the Central Limit Theorem)
"In distribution is the key!
(Alexander Matthews, NeurIPS 2019)
Matthews et al (2018),
generalisation of the result for multiple hidden layers:
- scaling priors correspondingly (normals with variance $\frac{1}{D}$, $\frac{1}{K_1}$, $\frac{1}{K_2}$ and so on)
As $K_1, K_2, ... , K_i$ all tend to infinity, the output layer $i$ is distributed as $N(0, C^{(i)}),$
where $C^{(i)}$ is recursively defined from $C^{(i-1)}.$
(Alexander Matthews, NeurIPS 2019)
Bayesian Decision Theory¶
Decision Theory¶
Given what we currently know about the world, how should we decide what to do, taking into account future events and observations that may change our conclusions?
Decision theory offers a way to automatically design and run experiments and to optimally construct clinical trials.
The exploration-exploitation trade-off¶
As we gather more information, we learn more about the world. However, the things we learn about depends on what actions we take.
For example, if we always take the same route to work, then we learn how much time this route takes on different days and times of the week. However, we don’t obtain information about the time other routes take. So, we potentially miss out on better choices than the one we follow usually. This phenomenon gives rise to the so-called exploration-exploitation trade-off.
The exploration-exploitation trade-off arises when data collection is interactive.
Bayesian Decision Theory¶
$\theta$ - parameter
$p(\theta)$ - prior on $\theta$
$A$ - all possible actions one can take
$l: A x \Theta \to R$ - loss function
The posterior risk of an action $a \in A$ is $$r(a) = E_\theta [ l(a, \theta) \mid D] = $$
i.e. the expectation of the loss
$$ = \int l(a, \theta) p(\theta \mid D) d \theta $$Loss functions:¶
- $l (\theta, \hat{\theta}) = (\theta - \hat{\theta})^2$. To compute the action, we need to minimize $$r(\hat{\theta}) = E [(\theta - \hat{\theta})^2 \mid D].$$ To derive the answer, take derivative of the integral with respect to $\hat{\theta}.$ And so $\hat{\theta} = E [\theta \mid D]$
$l (\theta, \hat{\theta}) = \mid \theta - \hat{\theta} \mid$ -> median of the posterior
$l (\theta, \hat{\theta}) = I_{\theta \ne \hat{\theta}}$ -> mode
Utility¶
Probability can be used to describe how likely an event is; utility can be used to describe how desirable it is.
Say again, frequentist or Bayesian?¶
Frequentist¶
- There are many ways to come up with a point estimator: the method of moments, etc. Our responsibility is to justify the choice of the estimator.
Bayesian¶
- There is no choice of the procedure of obtaining the estimate. Our choices are the prior, probability model, and (maybe) the loss function. Everything after that is computation!
How to chose a PPL?¶
- functionality (e.g. availability of discrete parameters)
- open to custom distributions (and samplers)
- preformance
- well documented
- abundant library of examples
- active (and supportive) community
An application from Pharma¶
an example of Julia + Turing being used for a publication
Take home messages¶
1. Advantages of Bayesian inference¶
- allows to use domain knowledge about the research question
- incorporate various sources of knowledge
- can answer a broad variety of question by computing desired qunatities from posterior samples
- can work with small data
2. Disadvantages of Bayesian inference¶
- computationally involved
- even running a model in a PPL can take a lot of time
- for MCMC moethods convergence is achieved, in theory, only at infinite number of iterations
3. PPLs¶
- modern approach to inference
- provide the iference "button"
- only the model needs to be formulated
- each PPL has its own syntaxis
Literuature and resources¶
learn Julia:
Bayesian Infernce:
- Robert, Casella, "Monte Carlo Statistical Methods"
- Robert, "The Metropolis–Hastings algorithm"
- Richard McEarlth, "Statistical ReThinking": https://xcelab.net/rm/statistical-rethinking/
- Stan examples and reference manual: https://github.com/stan-dev/example-models/wiki
- M. Betancourt, "Diagnosing biased inference with divergences": https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html
- M. Betancourt, "A Conceptual Introduction to Hamiltonian Monte Carlo"
- Model checking and diagnostics: https://pymc-devs.github.io/pymc/modelchecking.html
- Convergence diagnostics: https://rlhick.people.wm.edu/stories/bayesian_5.html
- Turing tutorials and documentation: https://turing.ml
- Bayesian Methods for Hackers: https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers
- "Predicting Drug-Induced Liver Injury with Bayesian Machine Learning", Williams, Lazic, Foster, Semenova, Morgan: https://pubs.acs.org/doi/pdf/10.1021/acs.chemrestox.9b00264
- Comparison of R/Stan, Python/PyMC3, Julia/Turing on the example of the Ordered Logostic Regression: https://medium.com/@liza_p_semenova/ordered-logistic-regression-and-probabilistic-programming-502d8235ad3f
using Pkg
Pkg.status()